load IntnlData MV;
MV(MV<=0) = NaN;
NumPort = 3; MinStocks = 50;

for cntry = [55 50 51]
    Line1 = cell(3,1); Line2 = cell(3,1);
    for xvar = 1:7
        X = XVARS(:,:,xvar);
        RPort = NaN(NumPort+1,NumPort+1,2,T);
        NPort = NaN(NumPort+1,NumPort+1,T);
        for t = TFirst-2:T-2
            R11_t = Return11(:,t);
            X_t1 = X(:,t+1);
            MV_t1 = log(MV(:,t+1));
            idx0 = isfinite(R11_t) & isfinite(X_t1) & isfinite(MV_t1) & CountryRegionDummy(:,cntry) & ~IsREIT & ~IsMicro(:,t,cntry);
            if sum(idx0) < MinStocks*NumPort^2 %number of stocks at this stage should be enough to populate all portfolios
                continue;
            end

            R11_t(idx0) = winsorize(R11_t(idx0),WinsorLevelX);
            X_t1(idx0) = winsorize(X_t1(idx0),WinsorLevelX);
            MV_t1(idx0) = winsorize(MV_t1(idx0),WinsorLevelX);
            RNext = Return(:,t+2);
            if WinsorizeYRet
                RNext(idx0) = winsorize(RNext(idx0),WinsorLevelY);
            end
            MVNext = MV(:,t+1);    MVNext(idx0) = winsorize(MVNext(idx0),WinsorLevelY);

            % neutralize country effects
            for cc = 1:49
                idx0_cc = CountryRegionDummy(:,cc) & idx0;
                if sum(idx0_cc) > 0
                    R11_t(idx0_cc) = R11_t(idx0_cc) - mean(R11_t(idx0_cc),'omitnan');
                    X_t1(idx0_cc)  = X_t1(idx0_cc) - mean(X_t1(idx0_cc),'omitnan');
                    MV_t1(idx0_cc)  = MV_t1(idx0_cc) - mean(MV_t1(idx0_cc),'omitnan');
                end
            end

            if xvar == 3 % convert to residual analyst
                y = X_t1(idx0); x = [ones(sum(idx0),1) MV_t1(idx0)];
                X_t1(idx0) = y - x*(x\y);
            end

            Xsort_t = sort(X_t1(idx0));
            NinPort = length(Xsort_t)/NumPort;
            BrkPts_m = Xsort_t([1 floor(NinPort*(1:NumPort))]); BrkPts_m(end) = Inf;
            if abs(BrkPts_m(1) - BrkPts_m(2)) < 1e-6 || abs(BrkPts_m(2) - BrkPts_m(3)) < 1e-6
                continue;
            end

            for m = 1:NumPort
                idx_m = X_t1>=BrkPts_m(m) & X_t1<BrkPts_m(m+1) & idx0;
                if sum(idx_m) < MinStocks*NumPort %number of stocks at this stage should be enough to populate all portfolios
                    continue;
                end

                R11sort_t = sort(R11_t(idx_m));
                NinPort = length(R11sort_t)/NumPort;
                BrkPts_n = R11sort_t([1 floor(NinPort*(1:NumPort))]); BrkPts_n(end) = Inf;
                for n = 1:NumPort
                    idx_nm = R11_t>=BrkPts_n(n) & R11_t<BrkPts_n(n+1) & idx_m;
                    if sum(idx_nm) >= MinStocks %ensure that number of stocks in this portfolio is enough
                        ret = RNext(idx_nm); mv = MVNext(idx_nm);
                        idx = isfinite(ret) & isfinite(mv);
                        ret = ret(idx); mv = mv(idx);
                        if sum(idx) >= MinStocks % ensure again that we have returns for enough stocks
                            RPort(n,m,1,t+2) = sum(ret.*mv)/sum(mv);
                            RPort(n,m,2,t+2) = mean(ret);
                            NPort(n,m,t+2) = sum(idx);
                        end
                    end
                end
            end
        end

        % take data only where corner portfolios are populated
        tmp = [squeeze(RPort(1,1,1,:)) squeeze(RPort(NumPort,1,1,:)) squeeze(RPort(1,NumPort,1,:)) squeeze(RPort(NumPort,NumPort,1,:))];
        idx = all(isfinite(tmp),2);
        RPort(:,:,:,~idx) = NaN; NPort(:,:,~idx) = NaN;

        for m = 1:NumPort
            RPort(NumPort+1,m,:,:) = RPort(NumPort,m,:,:) - RPort(1,m,:,:);
        end
        RPort(NumPort+1,NumPort+1,:,:) = RPort(NumPort+1,NumPort,:,:) - RPort(NumPort+1,1,:,:);
        for sample = 1:3
            switch sample
                case 1, T1 = TFirst; T2 = T;
                case 2, T1 = TFirst; T2 = THalf;
                case 3, T1 = THalf+1; T2 = T;
            end

            m = mean(RPort(:,:,:,T1:T2),4,'omitnan'); s = std(RPort(:,:,:,T1:T2),[],4,'omitnan'); n = sum(isfinite(RPort(:,:,:,T1:T2)),4); tstat = sqrt(n).*m./s; m = 1200*m;
            Line1{sample} = [Line1{sample} m(NumPort+1,NumPort+1,1)];
            Line2{sample} = [Line2{sample} tstat(NumPort+1,NumPort+1,1)];
        end
    end

    D = [Line1{1} NaN Line1{2} NaN Line1{3}; Line2{1} NaN Line2{2} NaN Line2{3}];
    if cntry == 55
        writematrix(D,FileName, 'Sheet', SheetName, 'Range','b17:x18', 'AutoFitWidth',false);
    elseif cntry == 50
        writematrix(D,FileName, 'Sheet', SheetName, 'Range','b20:x21', 'AutoFitWidth',false);
    elseif cntry == 51
        writematrix(D,FileName, 'Sheet', SheetName, 'Range','b23:x24', 'AutoFitWidth',false);
    end
end
